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Abstract: We report in this paper a fast and accurate algorithm for computing the exact spherical non- 
reflecting boundary condition (NRBC) for time-dependent Maxwell's equations. It is essentially based 
on a new formulation of the NRBC, which allows for the use of an analytic method for computing the 
involved inverse Laplace transform. This tool can be generically integrated with the interior solvers for 
challenging simulations of electromagnetic scattering problems. We provide some numerical examples 
Q ■ to show that the algorithm leads to very accurate results. 
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1. Introduction 

^ , The time-domain simulations, which are capable of capturing wide-band signals and modeling more 

p • general material inhomogeneities and nonlinearities, have attracted much attention lfT8l l9l[T2ll. A long- 

standing issue in many simulations resides in how to deal with the unbounded computational domain. 
Various approaches including the perfectly matched layer (PML) (cf. IS), the boundary integral methods 
3 ! (cf. 10), nonreflecting (absorbing or transpai^ent) boundary conditions (cf. ||8j|6]), and many others, have 

been proposed to surmount this obstacle. Although the use of exact nonreflecting boundary conditions 
(NRBCs) is desirable and beneficial, the practitioners are usually plagued with their complications and 
computational inefficiency. Indeed, these time-domain NRBCs are global in both time and space, and 
^ I involving Laplace inversion of special functions. It is worthwhile to highlight some works on efficient 

IjD • algorithms for exact NRBCs for acoustic wave equations, see e.g., |[T7l l2l lT3llT0l[T9l . However, there 

0> I has been significantly less study of the NRBCs for Maxwell's equations, where one finds the existing 

■ formulations (see e.g., lEJI^I) actually present a great challenge for evaluation. 

. I In this paper, we reformulate the NRBC for the three-dimensional Maxwell's equations, and extend 

■ the techniques for the NRBC of the acoustic wave equation in 121 [13 for computing it in a fast and 
accurate manner. It is important to point out that it is quite generic to integrate this sort of semi-analytic 
tool with any solver for the interior truncated problem (for example, the finite element/spectral element 
methods, and the boundary perturbation technique llT6l ). with the aid of the Spherepack 0] or certain 
hybrid mesh interpolation ifTTI . 

^ I Typically, we consider an electromagnetic scattering problem with a homogeneous background 

^ ■ transmission medium, and with a bounded scatterer D. Assume that the source current (or excita- 

tion source), other inputs and inhomogeneity of media are supported in a ball of radius b, that is, 
B := {jc € : < Z?). Then the analytic method of Laplace transform and separation of variables 
can be applied to solve the time-dependent Maxwell's system (exterior to B) with free source, homoge- 
neous initial data and the Silver-Muller radiation condition 

d,ET + cxxdtH = o{\xr^), t > 0; c = ll^[^, (1) 

where {E,H} are the electric and magnetic fields, x = x/\x\, and Ej = x x {E x x) is the tangential 
component of E. The electric permittivity e and magnetic permeability fi are positive constants. The 
underlying solution (cf. Hagstrom and Lau fl}) can be expressed in terms of vector spherical harmonic 
functions (VSHs) with the coefficients determined by the electric field on B: 

CO I 

^ = Z Z K^T^ + YT + ^Tr), ^ir = b, (2) 

1=1 in=-l 

where the VSHs |F™Jc, V5 y;", VsY'/' x x) are the orthogonal basis of {L\S))^ with S being the 

unit sphere (see e.g., fl4\ ). and {FJ") being the spherical harmonics as normalized in lITSl . Note that in 



Q, the exact NRBC is expressed as a system of E and H, which is actually equivalent to the formulation 
(cf. H) by using the VSH notation here: 



dtEr -cxxiV xE) = Th{E\ at r = b, (3) 
where the electric-to-magnetic (EtM) operator: 



''"['^1 = s Z Z * 4") yr + (T, , £«') Try (4) 



/=1 »!=-/ 

Here, pi and o"/ (termed as the nonreflecting boundary kernels (NRBKs)) are defined by 



Piit) = £- 



zkiiz) ^ ^ 



(t), (Ti{t) = £- 



k\{z) 

1 + Z + Z- 



'kliz) 



sb 

(0 with z = —, (5) 
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where £.~^ is the inverse Laplace transform (in ^-domain), and ki{z) is the modified spherical Bessel 
function defined by kj{z) = V27(^ii)^;+l/2(^)^ with Ki+1/2 being the modified Bessel function of the 
second kind of order 1+1/2 (cf. EOl ). The involved convolution is defined as usual: (/ * g){t) = 

£f{T)g{t-T)dT. 

Now, the central task is to compute the NRBC in Q-©. Notice that the electric field E at r = b 
is unknown as the NRBC sei^ves as the boundary condition for the interior problem. Here, we resort 
to the Spherepack |1| to communicate between the electric field and the VSH expansion coefficients. 
Thus, some hybrid mesh interpolation technique (cf. ifTTTl ') is necessary if the spatial discretization of 
the interior solver (e.g., the finite/spectral element methods) uses a different set of grids on the sphere. 
Thus, the critical issue becomes how to compute the NRBKs in (O, and temporal convolutions in ^ at 
any time t efficiently. This will be the topic of the following section. 

2. The Algorithm for Computing the NRBC 

The NRBK o"/ appears in the NRBC for the transient wave equation, which has an explicit formula 
(see ^ below) derived from the Residue theory (see e.g., ||7] [191). However, this analytic tool for 
inverse Laplace transform can not be applied to compute pi, since we lack information on the zeros of 
ki{z) + zk'jiz) (i.e., the poles of the integrand in the inverse Laplace transform), while that of ki(z) is 
available. In fact, there is no stable way to directly compute the NRBK p/. 

A. Alternative formulation of Tb[E]. 

Observe from dU that the EfM operator only involves the VSH expansion coefficients E^^] in 

Q. In fact, there holds the following relation between E^^^^ and E^'j^ : 



where \E'^J(S),e!"^^''^{s)] are Laplace transforms of {£'^^^^,£'1^1^^), respectively. The derivation of ^ is quite 
involved, so we will provide the proof in the extended paper. This leads to the following alternative 
formulation, from which the efficient algorithm stems. 

Theorem 1. The EtM operator Ti,[E] in @ can be reformulated as: 



00 / 



where cri is defined in Q, and oji is given by 

b , , _ , sb 



oJiit) = £' 



I b, , , sb 

z(l+z + z^)J(0 = -(^^;(0 + '^/(0)5(0), z = -, (8) 



with 5 being the Dirac delta function. 



B. Explicit formulas of the NRBKs cri and a*/. 

As already mentioned, the NRBK cr/ appears in the exact NRBC for the wave equation, and the 
explicit formula derived from the Residue theory (see e.g., ||7l[19l) is of the form: 

(^lit) = ' '^1' ^^0, (9) 

where are the zeros of ^/+i/2(z). Hence, it follows from ([8]) that 

/ / 
= J Yj-z]fei''^' + s{t) z'j, / > 1, t>0. (10) 

y=i 7=1 

Remark 1. We find from EOl that (see Figure [T] (left)): (i) Ki+1/2 has exactly / zeros, which appear in 
conjugate pairs and lie in the left-half of z-plane; and (ii) the zeros approximately sit along the boundary 
of an eye-shaped domain that intersects the imaginary axis approximately at ±i I, and the negative real 
axis at -la, where a w 0.66274. We point out that a practical algorithm in ifTOl could be used to find the 
zeros of A'/+i/2(z) for any / < 1000 accurately in negligible time. ■ 
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Figure 1: Left: distributions of the zeros of Ki+i/2iz)- Right: distributions of poles in ||2l(red color), and used 
zeros (blue color), which lies on the right of the vertical dashdot line -/3la with /3 = 0.4. 

Armed with the explicit formulas (l9l)-(fT0l). we can compute the NRBKs at any time. Observe that as 
the real part of Zj is negative, e*^-' becomes exponentially small when Zj is far away from the imaginary 
axis and t is slightly large. This motivates us to drop many insignificant zeros by using the algorithm 
described in |[T9l (see Figure [T](right)). 

C. Fast algorithms for temporal convolution. 

Observe from (l9]l-(fT0]i that the time variable t only appears in the exponentials. This allows for a fast 
recursive temporal convolution as shown in ||2]. More precisely, given a generic function g(t), we define 



and find 



fin r) ■= ei''' * g{t) - r ei'''''~''^g{T) dr, 
Jo 

Xt+At 



fit + At;r) = ei'-^'f{t;r) + 



(11) 



(12) 



where At is the time step size. We see that at each time step, the computation narrows down to computing 
the integral of the current intei^val [t,t + At]. This essentially eliminates the burden of history dependence 
induced by the temporal convolution. 

Using this notion, we deduce from (l9l)-(fTTI) that 



[T, * gm = ^ Z 4 e'^''fi-''>g{T)dT = ^ g z'j fit; z'j), 

[coi * gm = I Yj^z'jf \ ei''^^'-'\{T)dT + g{t) Yj^'j^l Yp'jfRt; z]) + g{t) ^ ^ 

■I Jo -I -I ;_ 1 



(13) 
(14) 



Given g at grids for time discretization, we only need to store {/(?; zyVj^i for previous steps to compute 
the convolutions at current time. It is optimal for storage requirement. 
D. Summary of the algorithm. 

For clarity, we summarize the whole algorithm and outline two strategies to further reduce the com- 
plexity related to large /. The algorithm is intended to compute 71, [£■](?) in ^ at t = t„ = nAt,n > 0, 
and on the colatitude-longitude grids adopted by the Spherepack HI from E on the same grids. 

Algorithm for computing Tt[E] in (|7]) 
Step 1. Use the Spherepack to compute {^J^,,, E^j^^} from E. 
Step 2. Compute the zeros {z'j} for / > 1. 

Step 3. Compute cr/ * Ef^j and a>/ * E^^^^ via ([T3])-(fT4]). 
Step 4. Use the Spherepack to compute ThlE] via ([7]). 



It is clear that the number of zeros to be used is determined by the truncation of the expansion ©. 
If I is large, we can adopt the strategies (i) dropping insignificant zeros (cf. |[T9l ) and (ii) compression 
algorithm (cf. Ill) to reduce the complexity in Step 3. Here, we outline the main idea. 

-z't I 

(i) . Since j becomes exponentially small when Zj is far away from the imaginary axis and t is 

slightly large, for some f^, > 0, we modify the NRBKs crj and ojj in (l9l)-(fT0l) as follows 

/ 

where if ^ {z^. : Re(z^.) > -pia) with a ~ 0.66742 and /? € (0, 1) (/? tunes the number of used 
zeros). We plot in Figure [T] (right) the used zeros for yS - 0.4. This can reduce the zeros from 
100 to 10 (see the marker 'x' on the right of vertical dashdot line) and leads to quite accurate 
approximation. We refer to the analysis in ||T9l on how to adjust /3 and ?g to achieve a good 
accuracy. 

(ii) . Alpert et al. [2J proposed a compression technique by a rational approximation of the NRBK in 

5-domain, which required to solve a nonlinear least square problem. This led to the approximate 
poles with d <s: I (see Figure [T] (right, marked by 'x') with a reduction from I = 100 to 

d = 12 and error tolerance 10^^). Correspondingly, o"/ and coi could be approximated by 

d d d 

cri{t)=^-Yj"''f ''''''' 'i'KO- ^2]c^-2j^'^'^' + <5(024 (16) 

;=i i=i ;=i 

where {cf^-ly^i are the coefficients occurring in the rational partial fraction: a'^jl{s - cz'j/b). Some 
samples of {a^,Zy}y^j are available from the website: http://faculty.smu.edu/thagstrom/sph6.txt. 

3. Numerical Results and Discussions 

In this section, we provide some numerical examples to show the accuracy of computing the NRBC. 
We also test a spectral-Galerkin with (second-order) Newmark's time integration for Maxwell's equa- 
tions in a spherical shell [a < \x\ < b}, where the NRBC is set at the outer spherical surface r = b. 

In the following tests, we generate the exact solution through the field: (x x E)\r=a = g, and with 
homogeneous initial data and source term. More precisely, we take 

^ ^ ^'( / ) 



where {gj'pgf2^ expansion coefficients in terms of {V5 Y"\ T'^}, of the the function 

/ 5 3 7 ■ 

g(6,(f>) = bi ^cos(;CjX2X3) + ly ,xl(^ sin (xixj) + l)" , a;i^0.5 cos (;c2a;3) + 1^" 



Note that x = (xi,X2,xi,) (with |jir| = 1) is the coixesponding Cartesian coordinates. Here, we compute 
{g"\,g'"2^ accurately by using the Spherepack lH]. 

We first test the accuracy of the algorithm for the NRBC dS}- Let be the truncated exact solution 
E (the included modes are 1 < / < Ne; \m\ < I). Define the error: 



e{b, t) - \\{dtE^ - c Jc X (V X £^)) ■ 



Ng,Ns 



(17) 



where 



denotes the discrete L^-norm associated with No x A^^ colatitude and longitude grids. 



Here, we take a - 2,c - 5,Ne - 40 and A^^ = INq. We aim to test the accuracy for computing T/;[£'^], 
so the differentiations in t and curl ai^e calculated analytically. In Table [T] we tabulate e{b, t) for different t 
and b (note: the magnitude of E is actually between 1 and 20, so the waves cross the artificial boundary). 
We see that in all cases, the computation of the NRBC is very accurate. 

Table 1 : Tfie error e{b, f) for different t and b 



t 


fe = 3 


b = 5 


b = 5 


b = 6 


1.0 


7.73246E-14 


2.67652E-14 


1.64591E-14 


1.10549E-14 


2.0 


7.01963E-14 


4.43072E-14 


4.30850E-14 


9.58167E-14 


4.0 


1.30672E-13 


7.07636E-14 


8.15754E-14 


9.56271E-14 


10.0 


3.35293E-13 


1.87063E-13 


2.32403E-13 


2.74720E-13 



Next, we set the NRBC as the boundary condition and solve the Maxwell's equation in curl-curl 
formulation with homogeneous initial conditions and free source in a spherical shell. In this case, we 
can expand the interior electric field in VSHs, and reduce the problem to a sequence of equations in 
radial direction. Then we solve the systems by using the spectral-Galerkin method in space and the 
second-order Newmark scheme in time. Moreover, we use the Richai^dson extrapolation to improve the 
time discretization to fourth-order. We refer to [|l9l for similar idea for the acoustic wave equations, and 
report the details in the extended version. 

Under the same setting of the reference solution and other data, we provide in Table |2J the discrete 
L^-norm errors (in space with sufficient resolution): Err(0 (Newmark scheme), Err^(0 (Richardson 
extrapolation), and the convergence order in time at different time t. As expected, we observe the second- 
order convergence for the Newmark scheme and the fourth-order convergence for the extrapolation. 



Table 2: Convergence of the Newmark scheme and Richardson extrapolation. 



t 


Af 


Err 


order 


Err"* 


order 




A/ 


Err 


order 


Err'* 


order 




5.00e-3 


2.3090E-2 




2.0719E-5 






5.00e-3 


1.6912E-2 




2.8491E-5 




0.5 


2.50e-3 


5.7684E-3 


2.001 


1.2870E-6 


4.009 


1.5 


2.50e-3 


4.2193E-3 


2.003 


1.7754E-6 


4.004 


1.25e-3 


1.4419E-3 


2.000 


8.0317E-8 


4.002 


1.25e-3 


1.0543E-3 


2.000 


1.1088E-7 


4.001 




l.OOe-3 


9.2276E-4 


2.000 


3.2912E-8 


3.998 




l.OOe-3 


6.7470E-4 


2.001 


4.541 lE-8 


4.009 




5.00e-3 


2.1041E-2 




2.6540E-5 






5.00e-3 


2.2117E-2 




2.3189E-5 




1.0 


2.50e-3 


5.2591E-3 


2.000 


1.6500E-6 


4.008 


2.0 


2.50e-3 


5.5142E-3 


2.004 


1.4467E-6 


4.003 


1.25e-3 


1.3147E-3 


2.000 


1.0299E-7 


4.001 


1.25e-3 


1.3776E-3 


2.001 


9.0376E-8 


4.001 




l.OOe-3 


8.4140E-4 


2.000 


4.2178E-8 


4.009 




l.OOe-3 


8.8159E-4 


2.000 


3.7016E-8 


4.000 



We have proposed in this paper an efficient algorithm for computing the spherical NRBC for the 
Maxwell's equations. This tool can be integrated well with various interior solvers in bounded domain 
for simulating scattering problems in many situations. We will report the works along this line in the 
forthcoming papers. 
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